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The late inspiral, merger, and ringdown of a black hole-neutron star (BHNS) system can provide 
information about the neutron-star equation of state (EOS). Candidate EOSs can be approximated 
by a parametrized piecewise-polytropic EOS above nuclear density, matched to a fixed low-density 
EOS; and we report results from a large set of BHNS inspiral simulations that systematically vary 
two parameters. To within the accuracy of the simulations, we find that, apart from the neutron- 
star mass, a single physical parameter A, describing its deformability, can be extracted from the 
late inspiral, merger, and ringdown waveform. This parameter is related to the radius, mass, and 
£ = 2 Love number, k2, of the neutron star by A = 2k2R^ /ZM^q, and it is the same parameter 
that determines the departure from point-particle dynamics during the early inspiral. Observations 
of gravitational waves from BHNS inspiral thus restrict the EOS to a surface of constant A in the 
parameter space, thickened by the measurement error. Using various configurations of a single 
Advanced LIGO detector, we find that neutron stars are distinguishable from black holes of the 
same mass and that A^''^ or equivalently R can be extracted to 10-40% accuracy from single events 
for mass ratios of Q = 2 and 3 at a distance of 100 Mpc, while with the proposed Einstein Telescope, 
EOS parameters can be extracted to accuracy an order of magnitude better. 

PACS numbers: 97.60. Jd, 26.60.Kp, 95.85.Sz 



I. INTRODUCTION 



Construction of the second-generation Advanced LIGO 
(aLIGO) detectors is underway, and will soon begin for 
Advanced VIRGO and LCGT, making it likely that grav- 
itational waveforms from compact binaries will be ob- 
served in this decade. Plans are also in development for 
the third generation Einstein Telescope (ET) detector 
with an order-of-magnitude increase in sensitivity over 
aLIGO. Population synthesis models predict that with 
a single aLIGO detector binary neutron star (BNS) sys- 
tems will be observed with a signal-to-noise ratio (SNR) 
of 8, at an event rate between 0.4 and 400 times per year 
and with a most likely value of 40 per year Black 
hole-neutron star (BHNS) systems are also expected, but 
with a more uncertain rate of between 0.2 and 300 events 
per year at the same SNR and with a most likely value 
of 10 events per year for a canonical 1.4 A^q-IO Mq sys- 
tem IJ. The expected mass ratios Q = Mbh/A^ns of 
BHNS systems are also highly uncertain and may range 
from just under 3 to more than 20 [21 [3]. 

A major goal of the gravitational-wave (GW) program 
is to extract from observed waveforms the physical char- 
acteristics of their sources and, in particular, to use the 
waveforms of inspraling and merging BNS and BHNS 
systems to constrain the uncertain EOS of neutron-star 
matter. During inspiral the tidal interaction between 
the two stars leads to a small drift in the phase of the 
gravitational waveform relative to a point particle sys- 
tem. Specifically the tidal field £ij of one star will in- 
duce a quadrupole moment Qij in the other star given 



by Qij = —^£ij where A"'^ is an EOS dependent quan- 
tity that describes how easily the star is distorted. A 
method for determining A for relativistic stars was found 
by Hinderer [1] ; its effect on the waveform was calculated 
to Newtonian order (with the relativistic value of A) by 
Flanagan and Hinderer [S] and to first post-Newtonian 
(PN) order by Vines, Flanagan, and Hinderer [6l[7]. This 
tidal description has also been extended to higher order 
multipoles [8, 9J. 

The detectability of EOS effects have been examined 
for both BNS and BHNS systems using this analytical 
description of the inspiral. For BNS systems, the de- 
tectability of A with aLIGO was examined for polytropic 
EOS [5j as well as a range of theoretical EOS commonly 
found in the NS literature for aLIGO and ET 'TU]. These 
studies considered only the waveform up to frequencies of 
400-500 Hz (~30-20 GW cycles before merger for 1.4 Mq 
equal-mass NSs). For this early part of inspiral, they 
find that the tidal deformability is detectable by aLIGO 
only for an unusually stiff EOS and for low neutron-star 
masses (< 1.2 Mq). ET on the other hand would have an 
order of magnitude improvement in estimating A, allow- 
ing ET to distinguish between different classes of EOS. 
For BHNS systems, using the recently calculated IPN 
corrections, Pannarale et al. [11] examined detectability 
for a range of mass ratios, finding that aLIGO will be 
able to distinguish between BHNS and binary black hole 



The tidal deformability for the Ith multipole is often defined in 
terms of the NS radius R and its dimensionless ^th Love number 
k( by Xi = f^2i—i)UG H^^'^^- Here we will discuss only the i = 2 
term so we write A := A2. 
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(BBH) systems only for low mass ratios and stiff EOS 
when considering the full inspiral waveform up to the 
point of tidal disruption. 

In sharp contrast to these analytic post-Newtonian 
results, analysis of just the last few orbits of BNS in- 
spiral from numerical simulations has shown that the 
NS radius may be extracted to a higher accuracy, of 
0(10%) [12], and this is confirmed by a study based on 
a set of longer and more accurate waveforms from two 
different codes [13]. In addition, comparisons between 
the analytical tidal description and BNS quasiequilib- 
rium sequences [T3] as well as long BNS numerical wave- 
forms [13 [H] suggest that corrections beyond the IPN 
quadrupole description are significant and substantially 
increase the tidal effect during the late inspiral. 

Numerical BHNS simulations have also been done to 
examine the dependence of the waveform on mass ratio, 
BH spin, NS mass, and the neutron-star EOS [T7H28] . 
However, an analysis of the detectability of EOS infor- 
mation with GW detectors using these simulations has 
not yet been done, and the present paper presents the 
first results of this kind. EOS information from tidal in- 
teractions is present in the inspiral waveform. For BHNS 
systems, however, the stronger signal is likely to arise 
from a sharp drop in the GW amplitude arising from 
tidal disruption prior to merger or, when there is negligi- 
ble disruption, from the cutoff frequency at merger [29] . 

We find from simulations of the last few orbits, merger, 
and ringdown of BHNS systems with varying EOS that, 
to within numerical accuracy, the EOS parameter ex- 
tracted from the waveform is the same tidal parameter 
A that determines the departure from point particle be- 
havior during inspiral; here A is a dimensionless version 
of the tidal parameter: 



Finally, we discuss future work in Sec. VII In the appen- 



A := GA 



GAh 



NS 



~ 3 ' \GMns, 



(1) 



where /c2 is the quadrupole Love number. 

The constraint on the EOS imposed by gravitational- 
wave observations of BHNS inspiral and merger is essen- 
tially a restriction of the space of EOS p = p{p) to a 
hypersurface of constant A, thickened by the uncertainty 
in the measurement (that is, a restriction to the set of 
EOS for which a spherical neutron star of the mass ob- 
served in the inspiral has tidal parameter A). We use 
a parametrized EOS based on piecewise polytropes [SH] . 
to delineate this region in the EOS space, but the result 
can be used to constrain any choice of parameters for the 
EOS space. 

In Sec. [n] we discuss the parametrized EOS used in 
the simulations. We give in Sec. [HI an overview of the 
numerical methods used and, in Sec. IVj a description of 
the waveforms from the simulations. We then discuss the 
analytical waveforms used for the early inspiral and is- 
sues related to joining the analytical and numerical wave- 
forms to create hybrids in Sec. IVj and we then estimate 



dices we discuss methods for numerically evaluating the 
Fisher matrix, and we provide instructions for generating 
effective one body (EOB) waveforms. In a second paper 
we will consider the detectability of EOS parameters for 
BHNS systems with spinning BHs. 

Conventions: Unless otherwise stated we set G = c = 
1. Base 10 and base e logarithms are denoted log and In 
respectively. We define the Fourier transform x{f) of a 
function x(t) by 



i(/)=/ x(t)e-2-/*di. 



and the inverse Fourier transform by 



(2) 



(3) 



II. PARAMETRIZED EOS 

To understand the dependence of the BHNS waveform 
on the EOS we systematically vary the free parameters of 
a parametrized EOS and then simulate a BHNS inspiral 
for each set of parameters. We choose the piecewise poly- 
tropic EOS introduced in Ref [30]. Within each density 
interval pi-i < p < pi, the pressure p is given in terms 
of the rest mass density p by 



Pip) = K,P' 



(4) 



where the adiabatic index is constant in each interval, 
and the pressure constant Ki is chosen so that the EOS 
is continuous at the boundaries pi between adjacent seg- 
ments of the EOS. The energy density e is found using 
the first law of thermodynamics, 



d- = —pd-. 
P P 



(5) 



the uncertainty in extracting EOS parameters in Sec. VI 



Ref. [nni uses a fixed low density EOS for the NS crust. 
The parametrized high density EOS is then joined onto 
the low density EOS at a density po that depends on the 
values of the high-density EOS parameters. The high- 
density EOS consists of a three-piece polytrope with fixed 
dividing densities pi = 10^^ '' g/cm^ and p2 = 10^^ g/cm^ 
between the three polytropes. The resulting EOS has 
four free parameters. The first parameter, the pressure 
Pi at the first dividing density pi , is closely related to the 
radius of a 1.4 Mq NS [31]. The other three parameters 
are the adiabatic indices {ri,r2,r3} for the three den- 
sity intervals. This parametrization accurately fits a wide 
range of theoretical EOS and reproduces the correspond- 
ing NS properties such as radius, moment of inertia, and 
maximum mass to a few percent |30j . 

Following previous work on BNS [H] and BHNS sim- 
ulations |17l I28j we use a simplified two-parameter ver- 
sion of the piecewise-polytrope parametrization and uni- 
formly vary each of these parameters. For our two pa- 
rameters we use the pressure pi as well as a single fixed 
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adiabatic index F = Fi = r2 = Fa for the core. The crust 
EOS is given by a single polytrope with the constants 
Kq = 3.5966 X 10^3 in cgs units and Fq = 1.3569 so that 
the pressure at 10^^ g/cra^ is 1.5689 x lO'^^ dyne/cm^. 
(For most values of pi, Fi, and F2, the central density of 
a 1.4 Mq star is below or just above p2, so the param- 
eter F3 is irrelevant anyway for BNS before merger and 
BHNS for all times.) 

We list in Table U the 21 EOS used in the simulations 
along with some of the NS properties. In addition, we 
plot the EOS as points in parameter space in Fig.[l]along 
with contours of constant radius, tidal deformability A, 
and maximum NS mass. The 1.93 Mq maximum mass 
contour corresponds to the recently observed pulsar with 
a mass of 1.97 ± 0.04 Mq measured using the Shapiro 
delay [32J . In this two-parameter cross section of the full 
four-parameter EOS space, parameters below this curve 
are ruled out. 
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the Einstein evolution equations in the BSSN formalism 
with the moving puncture gauge, and solves the hydro- 
dynamic equations with a high-resolution central scheme. 
The formulation, the gauge conditions, and the numerical 
scheme are the same as those described in Ref. [T7]. For 
the EOS, we decompose the pressure and energy density 
into cold and thermal parts as 



P = Pcoid +Pth , e = Ccoid + eth- 



(6) 



We calculate the cold parts of both variables using the 
piecewise polytropic EOS from p, and then the thermal 
part of the energy density is defined from e as eth = 
e — fcoid- Because eth vanishes in the absence of shock 
heating, eth is regarded as the finite temperature part. 
In our simulations, we adopt a F-law ideal gas EOS 



Pth = (Eth - l)eth, 



(7) 



to determine the thermal part of the pressure, and choose 
Fth equal to the adiabatic index in the crust region, Fq, 
for simplicity. 

In our numerical simulations, gravitational waves are 
extracted by calculating the outgoing part of the Weyl 
scalar ^'4 at finite coordinate radii ^ 4OOM0, and by 
integrating ^4 twice in time as 



h+{t)~ih^{t) 



t' 



dt' / dt"^4{t") 



(8) 



In this work, we perform this time integration with a 
"fixed frequency integration" method to eliminate un- 
physical drift components in the waveform |36| . In this 
method, we first perform a Fourier transformation of ^'4 
as 



FIG. 1: The 21 EOS used in the simulations are represented 
by blue points in the parameter space. For a NS of mass 
1.35 Mq, contours of constant radius are solid blue and con- 
tours of constant tidal deformability A are dashed red. Also 
shown are dotted contours of maximum NS mass. The shaded 
region does not allow a 1.35 Mq NS. 



III. NUMERICAL METHODS 

We employ BHNS binaries in quasiequilibruim states 
for initial conditions of our numerical simulations. We 
compute a quasiequilibrium state of the BHNS binary as 
a solution of the initial value problem of general relativ- 
ity, employing the piecewise polytopic EOS described in 
the previous section. The details of the formulation and 
numerical methods are described in Refs. [171 [33]. Com- 
putations of the quasiequilibrium states are performed 
using the spectral-method library LORENE |34| . 

Numerical simulations are performed using an 
adaptive-mesh refinement code SACRA ^35j. SACRA solves 



J —00 

Using this, Eq. ^ is rewritten as 



h+{t) - ihy,[t) 



) J 



(9) 



(10) 



We then replace 1//^ of the integrand with 1//q for 
I/I < /o, where fo is a free parameter in this method. 
By appropriately choosing /o, this procedure suppresses 
unphysical, low-frequency components of gravitational 
waves. As proposed in Ref. [35], we choose /o to be 
~ 0.8r7irio/27r, where is the initial orbital angular ve- 
locity and m(= 2) is the azimuthal quantum number. 



IV. DESCRIPTION OF WAVEFORMS 

Using the 21 EOS described in Table ^ we have per- 
formed 30 BHNS inspiral and merger simulations with 
different mass ratios Q — M^n/Mf^g and neutron star 
masses Mns- A complete list of these simulations is 
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TABLE I: Neutron star properties for the 21 EOS used in the simulations. The original EOS names [121 1171 128] are also listed. 
pi is given in units of dyne/cm^, maximum mass is in Mq, and neutron star radius R is in km. R, k2, and A are given for 
the two masses used: {1.20, 1.35} Mq. The values listed for logpi are rounded to three digits. The exact values used in the 
simulations can be found by adding log(c/cm s^^)^ - 20.95 ^ 0.00364 (e.g. 34.3 becomes 34.30364). 



EOS 


log Pi 


r 




-Rl.20 


^2,1-20 


Al.20 


-Rl.35 


^2,1.35 


Al.35 


p.3r2 


4 


Bss 


34.3 


2 


4 


1.566 


10.66 


0.0765 


401 


10.27 


0.0585 


142 


p.3r2 


7 


Bs 


34.3 


2 


7 


1.799 


10.88 


0.0910 


528 


10.74 


0.0751 


228 


p.3r3 





B 


34.3 


3 





2.002 


10.98 


0.1010 


614 


10.96 


0.0861 


288 


p.3r3 


3 




34.3 


3 


3 


2.181 


11.04 


0.1083 


677 


11.09 


0.0941 


334 


p.4r2 


4 


HBss 


34.4 


2 


4 


1.701 


11.74 


0.0886 


755 


11.45 


0.0723 


301 


p.4r2 


7 


HBs 


34.4 


2 


7 


1.925 


11.67 


0.1004 


828 


11.57 


0.0855 


375 


p.4r3 





HB 


34.4 


3 





2.122 


11.60 


0.1088 


872 


11.61 


0.0946 


422 


p.4r3 


3 




34.4 


3 


3 


2.294 


11.55 


0.1151 


903 


11.62 


0.1013 


454 


p.5r2 


4 




34.5 


2 


4 


1.848 


12.88 


0.1000 


1353 


12.64 


0.0850 


582 


p.5r2 


7 




34.5 


2 


7 


2.061 


12.49 


0.1096 


1271 


12.42 


0.0954 


598 


p.5r3 





H 


34.5 


3 





2.249 


12.25 


0.1165 


1225 


12.27 


0.1029 


607 


p.5r3 


3 




34.5 


3 


3 


2.413 


12.08 


0.1217 


1196 


12.17 


0.1085 


613 


p.6r2 


4 




34.6 


2 


4 


2.007 


14.08 


0.1108 


2340 


13.89 


0.0970 


1061 


p.6r2 


7 




34.6 


2 


7 


2.207 


13.35 


0.1184 


1920 


13.32 


0.1051 


932 


p.6r3 







34.6 


3 





2.383 


12.92 


0.1240 


1704 


12.97 


0.1110 


862 


p.6r3 


3 




34.6 


3 


3 


2.537 


12.63 


0.1282 


1575 


12.74 


0.1155 


819 


p.7r2 


4 




34.7 


2 


4 


2.180 


15.35 


0.1210 


3941 


15.20 


0.1083 


1860 


p.7r2 


7 




34.7 


2 


7 


2.362 


14.26 


0.1269 


2859 


14.25 


0.1144 


1423 


p.7r3 





1.5H 


34.7 


3 





2.525 


13.62 


0.1313 


2351 


13.69 


0.1189 


1211 


p.7r3 


3 




34.7 


3 


3 


2.669 


13.20 


0.1346 


2062 


13.32 


0.1223 


1087 


p.9r3 





2H 


34.9 


3 





2.834 


15.12 


0.1453 


4382 


15.22 


0.1342 


2324 



given in Table [TTj For the mass ratio Q = 2 and NS 
mass Mns = 1-35 Mq, we performed a simulation for 
each of the 21 EOS. In addition, we performed simula- 
tions of a smaller NS mass ((3 = 2, Mns = 1-20 Mq) 
and a larger mass ratio {Q = 3, A/ns — 1-35 Mq), 
in which we only varied the pressure pi over the range 
34.3 < log(pi/(dyne cm~^)) < 34.9 while holding the 
core adiabatic index fixed at F = 3.0. 

Two of the gravitational waveforms are shown in Fig. [2] 
below. The waveforms are compared with FOB BBH 
waveforms of the same mass ratio and NS mass which 
are also shown. Specifically we use the FOB formalism 
discussed in AppendixjB] The most significant differences 
begin just before the merger of the black hole and neu- 
tron star. For neutron stars with a small radius, the 
black hole does not significantly distort the neutron star 
which crosses the event horizon intact. As a result, the 
merger and ringdown of these waveforms are very simi- 
lar to the BBH waveform. However, a larger NS may be 
completely tidally disrupted just before merger resulting 
in a supressed merger and ringdown waveform. Disrup- 
tion suppresses the ringdown for two reasons related to 
the spreading of the matter: The ringdown is primarily 
a superposition of nonaxisymmetric quasinormal modes, 
dominated by the I = m = 2 mode, while the disrupted 



TABLE II: Data for the 30 BHNS simulations. NS mass is 
in units of Mq, and QqM is the angular velocity used in the 
initial data where M = Mbh + Mns • 



Q Mns 


EOS 


n 


oM 





Mns 


EOS QoM 


2 


1 


35 


p.3r2.4 





028 


2 


1.35 


p.6r3.3 0.025 


2 


1 


35 


p.3r2.7 





028 


2 


1.35 


p.7r2.4 0.025 


2 


1 


35 


p.3r3.0 





028 


2 


1.35 


p.7r2.7 0.025 


2 


1 


35 


p.3r3.3 





025 


2 


1.35 


p.7r3.0 0.028 


2 


1 


35 


p.4r2.4 





028 


2 


1.35 


p.7r3.3 0.025 


2 


1 


35 


p.4r2.7 





028 


2 


1.35 


p.9r3.0 0.025 


2 


1 


35 


p.4r3.0 





028 


2 


1.20 


p.3r3.0 0.028 


2 


1 


35 


p.4r3.3 





025 


2 


1.20 


p.4r3.0 0.028 


2 


1 


35 


p.5r2.4 





025 


2 


1.20 


p.5r3.0 0.028 


2 


1 


35 


p.5r2.7 





025 


2 


1.20 


p.9r3.0 0.022 


2 


1 


35 


p.srs.o 





028 


3 


1.35 


p.3r3.0 0.030 


2 


1 


35 


p.5r3.3 





025 


3 


1.35 


p.4r3.0 0.030 


2 


1 


35 


p.6r2.4 





025 


3 


1.35 


p.5r3.0 0.030 


2 


1 


35 


p.6r2.7 





025 


3 


1.35 


p.7r3.0 0.030 


2 


1 


35 


p.ers.o 





025 


3 


1.35 


p.9r3.0 0.028 
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matter is roughly axisymmetric as it accretes onto the 
black hole; and the accretion timescale of the spread-out 
matter is long compared to the periods of the dominant 
modes. 



The dependence of the waveform on the EOS can 
be seen more clearly by decomposing each waveform 
into amplitude A{t) and phase ^{t) with the relation 
h+{t) - ih-^{t) = ^(t)e-**(*). In Fig. [s) the amplitude 
as a function of time for each BHNS waveform is com- 
pared to a BBH waveform of the same value of Q and 
Mns- At early times, the waveform is almost identical 
to the BBH waveform. However, a few ms before the 
maximum amplitude is reached, the amplitude begins to 
depart from the BBH case. For each Q and Mns, this 
departure from the BBH waveform is monotonic in A. 
Neutron stars with large values of A merge earlier, and 
as a result the waveforms reach a smaller maximum am- 
plitude. The phase of each waveform is compared to that 
of the FOB BBH waveform $eob in Fig. |4] At early 
times the phase oscillates about the FOB phase due to 
initial eccentricity in the numerical waveform discussed 
in Sec. |VB[ At later times, closer to the merger, tidal in- 
teractions lead to a higher frequency orbit; this, together 
with correspondingly stronger gravitational wave emis- 
sion, means the BHNS phase accumulates faster than the 
FOB phase. This continues for 1-2 ms after the wave- 
form reaches its maximum amplitude (indicated by the 
dot on each curve). Eventually the amplitude drops sig- 
nificantly, and numerical errors dominate the phase. We 
truncate the curves when the amplitude drops below 0.01. 



The monotonic dependence of the waveform on A can 
again be seen in its Fourier transform h, shown in Figs.js] 
and [6j which is decomposed into amplitude and phase 
by h{f) = A(/)e-**(/). The predicted FOS dependent 
frequency cutoff in the waveform 29, is clearly shown in 
the amplitude^. Neutron stars that are more easily dis- 
rupted (larger A) result in an earlier and lower frequency 
drop in their waveform amplitude than NS with smaller 
A. The phase $(/) relative to the corresponding BBH 
waveform has a much smoother behavior than the phase 
of the time domain waveform. This feature will be use- 
ful in evaluating the Fisher matrix in Sec. |VI| The noise 
that is seen at frequencies above ~ 3000 Hz is the result 
of numerical errors in the simulation and has no effect on 
the error estimates below. 



^ Tidal disruption occurs after the onset of mass shedding of the 
neutron star. The frequency at the onset of mass shedding is 
usually much lower than that of tidal disruption for BHNS bi- 
naries |20| . In Ref. |29| . mass-shedding frequency was identified 
as the cutoff frequency but this underestimates the true cutoff 
frequency. See also Refs. [371 138| for a discussion of dynamical 
mass transfer. 



V. HYBRID WAVEFORM CONSTRUCTION 

Since our numerical simulations typically begin ^5 or- 
bits before merger, it is necessary to join the numerical 
waveforms to analytic waveforms representing the ear- 
lier inspiral. There is a substantial literature comparing 
analytic and phenomenological waveforms with numeri- 
cal waveforms extracted from simulations of BBH coa- 
lescence. For example, it has been shown that the 3.5 
post-Newtonian (TaylorT4) waveform agrees well with 
equal mass BBH waveforms up to the last orbit before 
merger [35.. For unequal mass systems, the FOB for- 
malism (see Ref. [40 for a review) has proven to be a 
powerful tool to generate analytic waveforms that agree 
with numerical simulations. Free parameters in the FOB 
formalism have been fit to numerical BBH waveforms to 
provide analytic (phenomenological) waveforms that ex- 
tend to the late, non-adiabatic inspiral as well as the 
ringdown. These FOB waveforms appear to be in good 
agreement with numerical BBH waveforms for mass ra- 
tios at least up to Q = 4 [¥!]• Although we have not ex- 
plored them in this context, other approaches have also 
been taken for constructing phenomenological inspiral- 
merger-ringdown waveforms j42H46| . 

For equal-mass BNS, Read et al. [12] compared the 
numerical BNS waveform during inspiral to a point par- 
ticle post-Newtonian waveform. Specifically, they used 
the 3.5 post-Newtonian (Taylor T4) waveforms matched 
on to the numerical waveforms to study the measurability 
of FOS parameters. They found that differences between 
the analytic and numerical waveforms become apparent 
4 — 8 cycles before the post-Newtonian coalescence time. 

The leading and post-l-Newtonian quadrupole tidal 
effects have recently been incorporated into the post- 
Newtonian formalism and used to compute corrections to 
the point-particle gravitational waveforms [SHZj- These 
post-Newtonian contributions along with a fit to the 
2PN tidal contribution have also been incorporated into 
the FOB formalism and compared to long simulations 
(^ 20 GW cycles), where they find agreement with the 
simulations to ±0.15 rad over the full simulation up to 
merger [16] . 

For the BHNS systems discussed here, we have 
matched the numerical waveforms to FOB waveforms 
that include inspiral, merger, and ringdown phases in- 
stead of post-Newtonian waveforms which are often not 
reliable during the last few cycles for higher mass ratios. 
This choice also allows us to use longer matching win- 
dows that average over numerical noise and the effects of 
eccentricity as shown in Sec. |V B| We have chosen to use 
the FOB formalism to generate inspiral-merger-ringdown 
waveforms, although we note that other phenomenologi- 
cal waveforms would probably work. For simplicity, and 
because it appears that an accurate description of the late 
inspiral dynamics just before merger requires 2PN tidal 
corrections |14H16j which are not yet known, we will use 
the FOB waveforms without tidal corrections. Our re- 
sults will therefore be lower limits on the measurability 




FIG. 2: h+ and \h\ = \h+ - ihx \ for BHNS waveforms for (Q = 2, Mns = 1-35 Mq) with two different EOS are represented 
by solid red and blue curves respectively. The softest EOS p.3r2.4 is on top and the stiffest EOS p.OFS.O is on bottom. An 
EOB BBH waveform (dashed) with the same values of Q and Mns is matched to each numerical waveform within the matching 
window Ti < t < Tf bounded by solid vertical lines. A hybrid EOB BBH-Numerical BHNS waveform is generated by splicing 
the waveforms together within a splicing window Si < t < Sf bounded by dotted vertical lines. The matching window is 12 ms 
long and ends at the numerical merger time tJv/' (time when the numerical waveform reaches its maximum amplitude), while 
the splicing window is 4 ms long and begins at the start of the matching window {Si —Ti). 



of EOS parameters since the EOS dependence is coming 
solely from the numerical waveforms. 



A. Matching procedure 

We use a method similar to that described by Read 
et al. [12] to join each of the numerical BHNS wave- 
forms to a reference EOB waveform, generating a hy- 
brid EOB-numerical waveform. Denote a complex nu- 
merical waveform by /iNR(i) = h^{t) — ihJ^(t) and an 
EOB waveform with the same Q and Mns by /ieob (^) = 
h^^{t) - ih^^^{t). A constant time-shift r and phase- 
shift $ can be applied to the EOB waveform to match 
it to a section of the numerical waveform by rewriting 
it as /iEOB(i — T)e^'*. We hold the numerical wave- 
form fixed because we must specify a matching window 
Tj < t < Tp, and as discussed below, there is only a 
small region of the numerical waveforms over which a 
valid match can be performed. Once the values of r 
and $ are determined, we will then choose to instead 



hold the EOB waveform fixed and shift the numerical 
waveform in the opposite direction by rewriting it as 
/i^hjft(t) = hys^it + T)e+**. This is done so that ah of 
the numerical waveforms with the same Q and Mns are 
aligned relative to a single fixed reference EOB waveform. 

Over a matching window Tj < t < Tp (bounded by 
solid vertical lines in Fig. [2|, the normalized match be- 
tween the waveforms is defined as 



m(r, $) 



Re [z{T)e"^] 
ctnro'eob(''' ) ' 



(11) 



where 



z{t)^ [ hNR{t)h*j,oB{t - t) dt (12) 

JTi 

and the normalizations for each waveform in the denom- 
enator are defined as 



4r= / \hMt)fdt (13) 
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FIG. 3: Amplitude of the complex waveform h — h+ — ihx ■ 
Dashed curves are EOB waveforms with the same Q and Mns • 
Matching and splicing conventions are those of Fig. [2] 



and 



c^eob(t) 



\hEOB{t - t)\'' dt. 



(14) 



Ti 



The time-shift r and phase $ are chosen to maximize 
the match m(T, for a fixed matching window. Ex- 
phcitly, the phase is determined analytically to be $ = 
'arg[z(r)]; plugging this result back into Eq. (11), the 



time-shift is given by maximizing |z(r)|/[(TNR0'EOB(''')] 
over r. As stated above, once r and $ are found we 
shift the numerical waveform in the opposite direction to 
generate hf^\t) = h^Rit + T)e+'*. 

A hybrid waveform is generated by smoothly turn- 
ing off the EOB waveform and smoothly turning on 



FIG. 4: Cumulative phase difference $ — $eob between 
BHNS waveform and EOB BBH waveform with the same Q 
and A^NS. The phase is defined by breaking up each com- 
plex waveform into amplitude and cumulative phase h+{t) — 



ihxit) = A{t)t 



The black point on each curve indi- 



cates the BHNS merger time tj^j defined as the time of max- 
imum amplitude The curve is truncated when the 



amplitude AD/M drops below 0.01. 
conventions are those of Fig. [2] 



Matching and splicing 



the shifted numerical waveform over a splicing window 
Si < t < Sf (bounded by dotted vertical lines in Fig. [2]) 
which can be chosen independently of the matching win- 
dow. As in Ref. [12j . we employ Hann windows 



Won{t) 
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FIG. 5: Weighted Fourier transform 2f^'^\h{f)\ of numer- 
ical waveforms where h = + hx)- Dot-dashed curves 
are EOB waveforms with the same Q and Mns . The left axis 
is scaled to a distance of 100 Mpc, and the noise sU^{f) 
for broadband aLIGO and ET-D are shown for compari- 
son. In each plot the numerical waveform monotonically ap- 
proaches the EOB waveform as the tidal parameter A de- 
creases. Matching and splicing conventions are those of Fig. [2] 
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FIG. 6: Cumulative phase difference $ — c&eob of the Fourier 
transform between BHNS waveform and EOB BBH wave- 
form of the same mass and mass ratio. The phase is defined 
by breaking up the Fourier transform h = \(h+ + hx) of 
each waveform into amplitude and cumulative phase h( f) = 
j4(/')e~'*(-'^'. Matching and splicing conventions are those of 
Fig. [2] 



The hybrid waveform is then viritten 

^hybrid (0 = 



/lEOB(t) t < Si 

Wos{t)hEOB{t) + Wo„(t)/l?Jf*(<) Si<t<SF 

h'^^'it) t > Sf 

(17) 

As shown in Fig. [2] we choose the start of the sphcing 
interval to be the same as the start of the matching win- 
dow Si = Ti and choose the end of the splicing window 
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to be Sp = Tj + A ms. It is also necessary to use these 
windows to smoothly turn on the hybrid waveform at low 
frequency when performing a discrete Fourier transform 
to avoid the Gibbs phenomenon. Unlike the case for BNS 
waveforms, it is not necessary to window the end of the 
hybrid waveform as the amplitude rapidly decays to zero 
anyway during the ringdown. 

For concreteness we define t = as the FOB BBH 
merger time tfP^ when the FOB waveform reaches its 
maximum amplitude. After matching to the FOB wave- 
form, the time when the numerical BHNS waveform 
reaches its maximum amplitude is t^p ■ 



B. Dependence on matching window 

Because the numerical BHNS waveforms are close but 
not identical to the FOB BBH waveform during the in- 
spiral and because there is some noise in the BHNS wave- 
forms, the time shift that maximizes the match depends 
on the choice of matching window. The matching window 
should exclude the first couple of cycles of the numeri- 
cal waveform during which time the simulation is settling 
down from the initial conditions. It should also exclude 
the merger/ringdown which are strongly dependent on 
the presence of matter. The window must also be wide 
enough to average over numerical noise and, as we shall 
see below, the effects of eccentricity in the simulations. 

The numerical merger time relative to the FOB 
BBH merger time tfp^ as a function of the end of the 
matching window Tp — t^f- provides a useful diagnos- 
tic of the matching procedure. Results for matching two 
Q = 2,Mns = I.35M0 waveforms with different equa- 
tions of state to an FOB waveform are shown in Fig. [7] 
The horizontal axis is the end time Tp of the matching 
window relative to the numerical merger time t^J^. For 
negative values, the matching window contains the BHNS 
inspiral only. For positive values, the matching window 
also contains part of the BHNS ringdown. The vertical 
axis is the location of the shifted numerical merger time 
t^/* after finding the best match. Four different window 
durations At = Tp — Tj are shown. The drift in the best 



fit merger time t^j^ most likely arises from the neglect of 
tidal effects in the FOB waveform which lead to an accu- 
mulating phase shift in the waveform, although it could 
also arise from numerical angular momentum loss from 
finite resolution of the simulations. Further work is in 
progress to understand this issue [13]. 

When the matching window duration is of order one 
orbital period or shorter, the time-shift oscillates as a 
function of Tp — t^p'. We attribute this effect to the 
eccentricity in the numerical waveform that results from 
initial data with no radial velocity. For larger matching- 
window durations, the effect of eccentricity is averaged 
out. 

To demonstrate concretely that the decaying oscilla- 
tions for At — 4 ms are the result of eccentricity, we 
matched an FOB BBH waveform with eccentricity to the 
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FIG. 7: Dependence of time shift between numerical and 
EOB waveform on the end time Tp — t^p' and width At of 
the matching window. Q — 2 and Mns = 1.35Mq for each 
waveform. The EOS used are p.3r2.4 (top panel), and p.QFS.O 
(bottom panel). The EOB waveform has zero eccentricity. 



equivalent zero eccentricity FOB BBH waveform. FOB 
waveforms can be generated with small eccentricity by 
starting the EOB equations of motion with quasicircular 
(zero radial velocity) initial conditions late in the inspi- 
ral. The result is shown in Fig. |8] for an FOB waveform 
with the same quasicircular initial conditions as the sim- 
ulation for the FOS p.3r2.4 shown in Fig. [T) The oscil- 
lations take exactly the form of those shown in Fig. [Tj 
except without the drift and offset. 

We estimate that the initial eccentricities in the sim- 
ulations used in this paper are Cq ~ 0.03. Decreasing 
the initial eccentricity by about an order of magnitude, 
possibly using an iterative method that adjusts the ini- 
tial radial velocity [47j, will remove this issue and allow 
one to determine the phase shift due to tidal interactions 
during the inspiral part of the simulation. 



VI. PARAMETER ESTIMATION 

The output of a gravitational- wave detector s{t) — 
n{t) -t- h{t) is the sum of detector noise n{t) and a possi- 
ble gravitational- wave signal h{t). Stationary, Gaussian 
noise is characterized by its power spectral density (PSD) 
5„(|/|) defined by 



{n{f)n*{n)^^5{f-f')SrS\f\) 



(18) 
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FIG. 8: Same as Fig.[7| but matching an eccentric EOB BBH 
waveform with the quasicircular initial condition A/fio = 
0.028 to a zero eccentricity EOB BBH waveform. 



The gravitational wave signal is given in terms of the two 
polarizations of the gravitational wave by 



h{t) = F+h+{t)+F^h^{t), 



(19) 



where -F+.x are the detector response functions and de- 
pend on the location of the binary and the polarization 
angle of the waves. We assume the binary is optimally 
located at the zenith of the detector and optimally ori- 
ented with its orbital plane parallel to that of the detec- 
tor. This condition is equivalent to averaging hj^. and hx 



Fx 



1/2). 



It is well known [JS] that the optimal statistic for de- 
tection of a known signal h{t) in additive Gaussian noise 
is 



(his) 



(20) 



where the inner product between two signals hi and h2 
is given by 



ihi\h2 



4Re 



Sn{f) 



df. 



(21) 



In searches for gravitational-wave signals from compact 
binary mergers, a parametrized signal h{t; 9"^) is known 
in advance of detection, and the parameters 6"^ must be 
estimated from the measured detector output s{t). The 
parameters 9^ of an inspiral are estimated by maximiz- 
ing the inner product of the signal s{t) over the template 
waveforms h{t;9^). In the high signal-to-noise limit, the 
statistical uncertainty in the estimated parameters 9^ 
arising from the instrumental noise can be estimated us- 
ing the Fisher matrix 



f dh 


dh \ 


\d9^ 


d9B ) 



(22) 



Note that 9 are the parameter values that maximize the 
signal-to-noise. The variance a\ = uaa = {{^S^)^) and 



covariance ctab = {A9^A9^) of the parameters are then 
given in terms of the Fisher matrix by 



{A9^A9^) ^ (F-i) 



AB 



(23) 



For hybrid waveforms, the partial derivatives in the 
Fisher matrix must be approximated with finite differ- 
ences. It is most robust to compute the derivatives of 
the Fourier transforms used in the inner product. We 
rewrite the Fourier transform of each waveform in terms 
of the ampl itude A and phase $ as exp[ln^ — i$] as given 



in Eq. (A9|. The derivatives d\nA/d9^ and 8^/89^ are 



then evaluated with finite differencing. More details of 
this and the other methods we tested are given in Ap- 
pendix |Aj 

In general, errors in the parameters 9^ are correlated 
with each other forming an error ellipsoid in parameter 
space determined by the Fisher matrix Tab- The un- 
correlated parameters that are best extracted from the 
signal are found by diagonalizing F^^. These new pa- 
rameters are linear combinations of the original parame- 
ters 9^. We focus attention below on the two parameters 
log(pi) and F, and fix all other parameters as follows. We 
use the masses and spins determined from the numeri- 
cal simulations and fix the time and phase shifts as de- 
termined during the hybrid waveform construction. We 
therefore construct the error ellipses in {log(pi),r} pa- 
rameter space and identify the parameter with the small- 
est statistical errors. We will leave an analysis of correla- 
tions due to uncertainty in masses and BH spin to future 
work. 



A. Broadband aLIGO and ET 

For the BHNS systems discussed here, the greatest de- 
parture from BBH behavior occurs for gravitational- wave 
frequencies in the range 500-5000 Hz. As a result, de- 
tector configurations optimized for detection of BHNS 
systems with low noise in the region below 500 Hz are 
not ideal for estimating EOS parameters. We therefore 
present results for the broadband aLIGO noise curve |35] 
and the ET-D noise curve [SOj shown in Fig. [g] The 
broadband aLIGO configuration uses zero-detuning of 
the signal recycling mirror and a high laser power, re- 
sulting in significantly lower noise above 500 Hz at the 
expense of slightly higher noise at lower frequencies. Sev- 
eral noise curves have been considered for the Einstein 
Telescope denoted ET-B [51], ET-C [52], and ET-D [50]. 
We will use the most recent ET-D configuration, and 
note that in the 500-5000 Hz range all of the ET config- 
urations have a similar sensitivity. The published noise 
curves, and those used in this paper, are for a single in- 
terferometer of 10 km with a 90° opening angle. The 
current ET proposal is to have three individual interfer- 
ometers each with a 60° opening angle. This will shift 
the noise curve down appoximatcly 20% jSOj . 

In Figs. [T0| and [TT] we show the resulting 1-a error el- 
lipses in the 2-dimensional parameter space {log(pi),r} 
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FIG. 9: Noise PSD for broadband aLIGO (dashed blue), 
ET-D (dotted red) and various configurations of narrowband 
aLIGO (solid black). The minima of the narrowband config- 
uration are labeled /r 



for an optimally oriented BHNS with Q — 2 and Mns = 
1.35M0 at a distance of 100 Mpc. Surfaces of constant 
A^/^ and NS radius, which are almost parallel to each 
other, are also shown. One can see that the error ellipses 
are aligned with these surfaces. This indicates that, as 
expected, A^/^ is the parameter that is best extracted 
from BHNS gravitational-wave observations. Because 
Ai/5 and R are so closely aligned we will use these two 
parameters interchangeably. 
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FIG. 10: Two 1-(T error ellipses for broadband aLIGO. 
Evenly spaced contours of constant A^^^ are also shown. Each 
ellipse is centered on the estimated parameter denoted by 
a X. The semimajor axes are significantly longer than the 
width of the figure, so each ellipse appears as a pair of par- 
allel lines. Matching and splicing conventions are those of 
Fig-i 

As mentioned above, there is some freedom in con- 
struction of the hybrid waveforms. The size and orien- 




FIG. 11: l-o" error ellipses for ET-D. Evenly spaced con- 
tours of constant A^''^ (i?) are also shown on top (bottom). 
Matching and splicing conventions are those of Fig. [2] 



tation of the error ellipses also depend on the details of 
this construction. We find that as long as the matching 
window is longer than approximately four gravitational- 
wave cycles to average out the effects of eccentricity and 
does not include the first two gravitational-wave cycles, 
the orientation of the error ellipses does not change sig- 
nificantly. As expected, the size of the ellipses decreases 
as more of the numerical waveform is incorporated into 
the hybrid waveform. We therefore adopt the last 12 ms 
before merger of each numerical waveform as the match- 
ing window and the first 4 ms of the matching window 
for splicing as shown in Fig. [2] 

We have emphasized that, to within present numer- 
ical accuracy, the late-inspiral waveform is determined 
by the single parameter h^l^ . This implies that, by us- 
ing countours of constant A in the EOS space, one could 
have obtained the constraint on the EOS, summarized in 
Figs . [To] and pT] by varying only a single EOS parameter. 
For the simulations with other mass ratios and neutron 
star masses, we have used as our single parameter log(pi) 
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and not T because countours of constant pi more closely 
coincide with contours of constant A and because A is 
a one to one function of log(pi) throughout the param- 
eter space. The one-parameter Fisher matrix can then 
be evaluated with finite differencing using the waveforms 
and values of A at two points in EOS parameter space 
with different log(pi). 
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The uncertainties in A^/^ and R are shown in Figs 
and [13] for broadband aLIGO and for ET respectively. 
The uncertainty in these quantities is ~ 10-40% for 
broadband aLIGO and ~ 1-4% for ET-D. The uncer- 
tainties for the higher mass ratio Q = 3 are somewhat 
larger than for Q — 2, but not significantly so. It is not 
clear how rapidly the uncertainty in A^/^ and R will in- 
crease as the mass ratio is increased toward more realistic 
values. On the one hand the tidal distortion is likely to 
be much smaller for larger Q. On the other hand the 
overall signal will be louder, and the merger and ring- 
down will occur at lower frequencies where the noise is 
lower. Additional simulations for higher Q are needed to 
address this question. 



2.0r 



0.20 r 



0.15 



<0.10- 

b 



0.05 

o.o^E 

0.5p 
0.4- 
-0.3- 



b 0.2- 



0.1 - 



0.0 1 



10 



aQ=2, MNs=1.20Mg 
□ Q=2, MNs=1.35Mg 
*Q=3, Mns=1.35M 



3.0 



3.5 



4.0 



4.5 



5.0 



.1/5 



aQ=2, /MNs=1.20/Wg 
□ Q=2, /M|,js=1.35/Mg 
*Q=3, /Mns=1.35/W„ 



□ 
□ 



□ 
▲ 



* 

□ 



11 



12 



13 
R (km) 



14 



15 



16 



1.5- 



aQ=2, /MNs=1.20/Mg 
□ Q=2, /Wns=1.35/W^ 
* Q=3, /Wns=1 .35M 



FIG. 13: Same as Fig. [12] but with the ET-D noise PSD. 
Error estimates are an order of magnitude smaller than for 
broadband aLIGO. 
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FIG. 12: 1-a uncertainty cr^i/s and cr^j as a function of the 
parameters A^''^ or R for the broadband aLIGO noise PSD. 
Matching and splicing conventions are those of Fig. [2] 



B. Narrowband aLIGO 

The presence of a signal-recycling cavity in the aLIGO 
instruments will allow them to be tuned to have improved 
narrowband sensitivity at the expense of bandwith. Two 
parameters control the narrowband capabilities of the in- 
struments [53-55j: the signal recycling mirror transmis- 
sivity effectively sets the frequency bandwidth of the in- 
strument, while the length of the signal recycling cavity 
(or equivalently the signal-recycling cavity tuning phase) 
controls the central frequency fn of the best sensitivity. 
By tuning one or more of the aLIGO detectors to oper- 
ate in narrowband mode, it may be possible to improve 
estimates of the EOS parameters. 

We have examined several narrowband tunings with 
central frequencies that vary between approximately 
fj^ = 500 Hz and 4000 Hz. These noise curves use a sig- 
nal recycling mirror transmissivity of 0.011 and a signal- 
recycling cavity tuning phase ranging from 10° down to 
1°, and were generated using the program gwinc [56]. 
Three of these noise curves are shown in Fig.|9] In Fig. [14] 
we plot the 1-a uncertainty in NS radius aji as a function 
of the narrowband central frequency /r. For the wave- 
forms considered in this paper the optimal narrowband 
frequency is in the range 1000 Hz < //j < 2500 Hz and 
depends on the mass ratio, NS mass, and EOS. Narrow- 
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band configurations usually give smaller errors than the 
broadband configuration if fn happens to be tuned to 
within a few hundred Hz of the minimum for that BHNS 
event. In Ref. ^S7], Hughes discussed a method for de- 
termining the best frequency fn to tune a narrowbanded 
detector to extract an EOS dependent cutoff frequency 
from a sequence of identical BNS inspirals. While this 
technique is not directly applicable to BHNS systems, 
which have different masses and spins, a similar approach 
could be used to combine multiple BHNS observations. 
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FIG. 14: 1-a uncertainty in R for different configurations of 
narrowband aLIGO and for different EOS. fn defines the fre- 
quency where Sn is a minimum as shown in Fig. [9] Horizontal 
lines on the left and right indicate the corresponding 1-a er- 
rors for broadband aLIGO and ET-D respectively. Matching 
and splicing conventions are those of Fig. [2] 



VII. DISCUSSION 



Results 



Using a large set of simulations incorporating a two- 
parameter EOS, we have found that the tidal deforma- 
bility A^/^, or equivalently the NS radius i?, is the pa- 
rameter that will be best extracted from BHNS wave- 
forms. These parameters can be estimated to 10-40% 
with broadband aLIGO for an optimally oriented BHNS 
binary at 100 Mpc. The narrowband aLIGO configu- 
ration can do slightly better if it is tuned to within a 
few hundred Hz of the ideal frequency for a given BHNS 
event. The proposed Einstein Telescope will have an 
order-of-magnitude better sensitivity to the EOS param- 
eters. 

Although we have used a particular EOS parametriza- 
tion to show that A is the parameter that is observed 
during BHNS coalescence, this result can be used to con- 
strain any EOS model — an EOS based on fundamental 
nuclear theory in addition to a parametrized phenomeno- 
logical EOS. In particular, several parametrizations have 
recently been developed, including a spectral represen- 
tation |58j . a reparametrization of the piecewise poly- 
trope [53] , and a generalization that also includes nuclear 
parameters [50] . 

The results presented here can be compared with re- 
cent work to determine the mass and radius of individual 
NS in Type-1 X-ray bursts. Ozel et al. [HI] have obtained 
mass and radius measurements from several systems by 
simultaneously measuring the flux F, which is likely close 
to the Eddington value, and the blackbody temperature 
T during X-ray bursts of systems with accurately deter- 
mined distances. During the burst, the emission area 
of the photosphere F/{aT'^) expands, contracts, then 
reaches a constant value, and Ozel et al. have argued 
that the final area corresponds to that of the NS surface. 
They obtain estimates of NS mass and radii with 0(10%) 
1-a uncertainty. Steiner et al. |60j have also considered 
these systems, but argue that the final emission area does 
not necessarily correspond to that of the NS surface, and 
as a result obtain slightly smaller NS radii and larger un- 
certainties in the mass and radius. These radius error 
estimates are slightly smaller than those for the BHNS 
systems we have considered at 100 Mpc. However, we 
note that binary inspiral observations are subject to less 
systematic uncertainty due to questions of composition of 
the photosphere and associating it with the NS surface. 

The uncertainty in NS radius for the merger and ring- 
down of BHNS systems examined here is of roughly 
the same size as that found for the last few orbits up 
to merger of BNS systems at the same 100 Mpc dis- 
tance [m [T3]. BNS inspirals, however, will likely oc- 
cur more frequently, and, including a tidally corrected 
inspiral-numerical hybrid, BNS systems are likely to have 
uncertainties that are smaller than BHNS systems by a 
factor of a few. Considering the post-merger phase for 
BNS waveforms may also provide additional information. 
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Expected NS masses in both BNS and BHNS systems are 
slightly smaller than those measured for X-ray bursters 
which have accreted matter from their companion, so 
BNS and BHNS GW observations may complement X- 
ray burst observations by better constraining the lower 
density range of the EOS which is not well constained 
from X-ray burst observations jSHl ISI] • 

B. Future work 

We have used in this paper several simplifications and 
conventions which can significantly effect the accuracy to 
which EOS parameters can be extracted. We list them 
below and describe how changing them would effect the 
parameter error estimates. 

1. Finite length of numerical waveforms 

The BHNS waveforms used here include only the 
last ~10 GW cycles of inspiral as well as the merger 
and ringdown, of which the first few cycles of in- 
spiral are unreliable due to inexact initial data. 
Matching the numerical waveform to a tidally cor- 
rected inspiral waveform instead of just the point- 
particle waveform will increase the overall depar- 
ture from point-particle behavior by (i) creating a 
phase shift during the early inspiral, and more im- 
portantly (ii) adding to the phase of the late inspi- 
ral and merger the accumulated phase shift from 
the early inspiral ~ a phase shift that is not already 
included in the stronger signal of the late inspiral. 
The tidal corrections are now known up to IPN or- 
der. For simulations with the current number of 
orbits, however, it appears that higher order tidal 
corrections will be needed to fully describe the late 
inspiral where matching to numerical waveforms is 
done. We leave the issue of generating tidally cor- 
rected inspiral-numerical hybrid waveforms to fu- 
ture work. 

2. Event rates 

Estimates of the detectability of EOS parameters 
in BNS systems are often given for an event at a 
distance of 100 Mpc, and we have used the same 
convention here to state the results above. The 
relevant event rate is, therefore, the expected num- 
ber of detected events that will have an effective 
distance Dcs < 100 Mpc. (The effective distance 
DeS depends on the location of the binary and 
its inclination relative to the detector. For an 
optimally oriented and located binary, one finds 
D = D^ff while typically D < D^g.) The aLIGO 
inspiral rates for BNS systems are highly uncertain 
with {low, most likely, high} estimates of {0.01, 1, 
10} Mpc"3 Myr-i [T] or {0.004, 0.4, 4} yr'^ with 
effective distance Dcs < 100 Mpc. Rates are even 
more uncertain for BHNS systems with rate esti- 
mates of {0.0002, 0.01, 0.4} yr-i with effective dis- 
tance Dcs < 100 Mpc [T]. Since the uncertainty 



in EOS parameters scales linearly with distance 
[cTyYi/s = (7y\^i/5 iooMpc(-C'/100Mpc)] and the event 
rate scales as D^, the estimated detection rates of 
systems with effective distance DcB < 400 Mpc are 
{0.01, 1, 30} yr^^ with a four-fold increase in un- 
certainty of A^/^. Fortunately, for iVobs identical 
events and N^ct identical detectors, the uncertainty 
also scales as ctavs /V^obs-^dct- 

3. Expected NS masses and mass ratios 

The simulations we used included realistic mass 
neutron stars of 1.2 and 1.35 Mq. On the other 
hand, black hole masses are expected to be many 
times larger [5], with likely mass ratios closer to 
Q ~ 7 (for the canonical 10 M0-1.4 Mq system) 
than the Q — 2 and 3 systems we examined here. 
Additional simulations for mass ratios of 4 and 5 
are in progress. 

4. Spinning BH 

In this paper we have not examined the effect of a 
spinning BH. The analytic results of Ref. [11 indi- 
cated that spin does not significantly improve the 
sensitivity to A for the inspiral up to the point 
of tidal disruption. However, numerical simula- 
tions [24l |27l [28] have shown that spin can strongly 
affect the dynamics near tidal disruption and the 
amount of matter left over in an accretion disk. We 
have performed several tens of simulations of non- 
precessing BHNS systems with spinning BH with 
various BH spins, mass ratios, NS masses, and EOS 
parameters, and an analysis of how BH spin affects 
the detectability of EOS parameters will be the sub- 
ject of the next paper. 

5. Correlations between parameters 

In our Fisher analysis we have assumed that the 
mass ratio, NS mass, and BH spin will be deter- 
mined to sufficient accuracy during the inspiral to 
separate them from EOS effects during the merger 
and ringdown. A full Fisher analysis using all of the 
BHNS parameters should be done to find the extent 
to which uncertainties in the other parameters alter 
measurability estimates of the EOS parameters. 

Because BHNS waveforms smoothly deviate from cor- 
responding BBH waveforms as A increases, it is likely 
that one can find a good analytical approximation for the 
full inspiral, merger, and ringdown waveform by modify- 
ing analytical BBH waveforms. Accurate waveforms for 
non-spinning BBH systems using the FOB approach have 
been developed IH] [SS] and work to find FOB waveforms 
for spinning BBH systems is in progress |64[ I65j . Tidal 
interactions have also been incorporated into the FOB 
approach for BNS systems with good agreement with 
the inspiral waveform from numerical simulations when 
parametrized 2PN tidal interactions are fit to the numer- 
ical waveform |151 116j . Another approach is to use phe- 
nomenological waveforms that fit the frequency domain 
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post-Newtonian inspiral waveform to a phenomenological 
merger and ringdown for both spinning and non-spinning 
BBH systems [32]. Both of these approaches may work 
for generating full analytic BHNS waveforms as well. A 
complete description of the BHNS waveform will likely 
include corrections for the I — 3 tidal field and other 
higher order corrections. However, it is not clear given 
the current set of simulations that these effects would 
be observable with either aLIGO or a third generation 
detector such as ET. 
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Appendix A: Numerically evaluating the Fisher 
matrix 

When an analytical representation of a waveform is 
not available, the partial derivatives in the Fisher matrix 



Eq. ( 22 1 must be evaluated numerically. There are sev- 



eral possible methods one can use, and we will examine 
their accuracy below. 



1. Finite differencing of h{t;9) 

The simplest method, and that used in Ref. [12] . 
is straightforward finite differencing of the signal h = 
F-^-h^ + Fxhx ■ For example, for five waveforms with val- 
ues of an EOS parameter 9 given by {9^2,6-1,(^0, 61,62} 
with equal spacing A9, the three and five point central 
differences are given by 



dh 
d9 

A2h 

A6 
dh 

d9 

Aih 

~d9~ 



A.h 



Ad 

-lh{t;9^i) 



0{{A9f), where 
\h{t-9i) 



AaH 



A9 

hKt;6 



A9 

0{{A9f), 



where 



\h{t-9_ 



(Al) 



lh{t■9^)~ ^h{t-92) 



A9 



(A2) 



This finite differencing method is useful when waveforms 
differ only slightly: at each time t, on the scale A6 the 
function h{t] 9) is well approximated by the low order 
interpolating polynomials used to generate the finite dif- 
ferencing formulas. 
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This assumption fails when the waveforms used in the 
finite differencing are significantly out of phase with each 
other'^. The tidal interaction leads to a monotonically 
accumulating phase difference relative to a BBH wave- 
form, implying that at a fixed time t the function h{0;t) 
is an oscillating function of 0. Now if an oscillating func- 
tion h[^{9)] = cos[$(0)] has wavenumber k = ^'{9) that 
varies slowly compared to $, then h'{9) is better approxi- 
mated by - sin($)A$/A6' than by A cos[$(6l)]/A6l. The 
assumption that k is slowly varying is k' <^ k^, k" <C k^ ^ 
and the error in, for example, each of the two second- 
order discretizations is given to order A0^ by 



dh 
d0 



dh _ 
d9 ~ 

sin[$(6i) 



= h{9)[lk^ + Oikk',k")]Ae^, 



A9 



(A3) 



with the error in the second expression much smaller than 
that in the first. We consider two ways to take advantage 
of this difference in accuracy. 



2. Finite differencing of amplitude and phase 

The first is to decompose each complex waveform into 
an amplitude A and accumulated phase $ 



h+it; 9) - ihy, {t- 9) = A(t; 0)e-'*(*'^\ 



(A4) 



where the accumulated phase of each waveform is a con- 



tinuous function defined by $ 



arg(/i_, 



for some integer n, and at the starting time ti the ac- 
cumulated phase of each waveform is chosen to be on 
the branch n — Q. The advantage of this method is 
that, at a fixed time, the functions A{t] 9) and 9) are 
non-oscillatory functions of 9 even when the accumulated 
phase difference between two waveforms is significantly 
more than a radian. 

With this decomposition the gravitiational wave signal 

is 

h{t]9) = A{t-9){F+coii^{t;9) + F^ sin$(t;6')), (A5) 
and the derivative of h is approximated by 
dh AA, 



d9 A9 



(F+ cos $ + -Fx sin <i>) 



-A{-F+ sin $ + Fx cos <&) 



A$ 
A0' 



(A6) 



If an intermediate waveform is not available to provide 
the functions A{t; 9o) and 9o), they can be evaluated 
by e.g. A{t; 9o) = {A{t; + A{t; 9,))/2. 



We have found that this method works reasonably well 
for the inspiral waveform. If, however, the amplitude of 
one of the numerical waveforms drops to zero, then the 
phase of the waveform becomes undefined. Because the 
amplitude of the numerical BHNS waveforms fall to zero 
at different times for different EOS, as shown in Fig. [3) 
the derivative d^/d9 becomes meaningless towards the 
end when the average amplitude is still nonzero. It is 
likely one could work around this difficulty. However, we 
choose instead to use another more robust method. 



3. Finite differencing of Fourier transform 

Because we will need to calculate the Fourier trans- 
form of the derivative dh/d9 to find the Fisher matrix, 
we first Fourier transform each waveform and then eval- 
uate the numerical derivative. Since the derivative d/d9 
commutes with the Fourier transform, the Fisher matrix 
can be written explicitly as 



/ dh 
\d9^ 



dh 
d9B 



= 4Re 



fl dh dhZ 
09^ de^ 



fi 



Snif) 



df, (A7) 



where the contribution to the integral below fi and above 
// is negligable. 

As in the second method we break up each Fourier 
transformed waveform into amplitude A{f;9) and accu- 
mulated phase $(/;6') 



Hf;d)=A{f;9)e-^-^(f--'\ 



(A8) 



where the phase of each waveform at fi is on the n — 
branch cut. As demonstrated by Figs. [5]and|6j both the 
amplitude and phase are non-oscillatory functions of 9 
at a fixed frequency /, and can be well approximated by 
a low-order polynomial. In contrast to the accumulated 
phase of the complex numerical waveform h^ — ihx , the 
accumulated phase of the Fourier transform of the strain 
h is always well defined for numerical BHNS waveforms 
in the frequency range fi to //. 

Finally, we find that one obtains better accuracy by 
differentiating In A instead of A, decomposing h as 



,lnA(/;9)-i$(/;e) 



(A9) 



The derivative is now approximated by 

89 \ A9 A9 J ' 

interpolating when needed to evaluate In A and $ at the 
midpoint. 



(AlO) 



The dephasing of numerical waveforms is even more significant 
for BNS inspiral. We believe that Ref. |12j which used this 
method underestimated the derivatives in some cases by a factor 
of ~2 or more, and thus overestimated the uncertainty in EOS 
parameters by the same factor. 



4. Parameter spacing 

Finally, we note that the EOS parameter spacing must 
be carefully chosen. If two waveforms are too close in 
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ds^ = -A{r)dt^ + B{r)dr^ + r^{dB 



parameter space, the error in each waveform will domi- 
nate over the truncation error due to finite differencing. 
The most significant source of this error comes from the 
spurius oscillations in the amplitude of the Fourier trans- 
form in the frequency range ^ 500-800 Hz (see Fig. [5| 
that result from joining the FOB and numerical wave- 
forms which are not exactly the same in the matching 
window. We find that the integrand of the Fisher ma- 
trix is often erratic in the range ~ 500-800 Hz when 
using the smallest parameter spacing available. How- 
ever, when the spacing is increased, the integrand is 
smoother in this frequency range and its contribution 
to the integral is significantly reduced. For the mass ra- 
tio (5 = 2, we find that a spacing between waveforms of ^ 
Alog(pi/(dyne cm'^)) = 0.1 for the first EOS parame- ^(") = ^5 [l-2u+2i 
ter is often sufficiently large to reduce this problem, while 
a spacing of AF — 0.6 for the second EOS parameter is 
the minimum spacing one can use. For Q = 3, we have 
found that a spacing of A log(pi/(dyne cm~^)) > 0.2 is 
necessary to reduce this problem. 

In addition, if the EOS parameters of two waveforms 
lie near the same degenerate contour where waveforms 
are identical (e.g. a contour of constant A which can be 
nearly identical to a line of constant F) , the error in each 
waveform will again dominate the truncation error even 
if the EOS parameters are widely spaced. For our two- 
dimensional EOS parameter space, this problem can be 
solved by transforming the parameter space such that 
points that originally formed a x pattern now form a 
-|- pattern, and in the transformed parameter space the 
new axes are not alligned with a degenerate contour. The 
Fisher matrix can be calculated in the transformed pa- 
rameter space then transformed back to the original pa- 
rameter space. 

We find that as long as these two requirements are met, 
the uncertainties in CT;^i/5 and have only an 0(20%) 
fractional dependence on the EOS parameter spacing. 
However, the dependence of the orientation of the error 
ellipses on the EOS parameter spacing does not allow one 
to distinguish between A and R as the best extracted pa- 
rameter. 



1. Hamiltonian dynamics 

In the EOB formalism the two-body dynamics are re- 
placed by a test particle of reduced mass /i — 77117712/ M 
moving in a modified Schwarzschild metric of total mass 
M = mi + m2 given by^ 



2 ' sin2 



'■)■ (Bl) 



The metric potentials A and B can be calculated from 
post-Newtonian theory. The first function is 

94 4l7r2\ 4 5 g 

— — 1 lyu +azvu +aevu 

(B2) 

where u = 1/r, ly = fi/M is the symmetric mass ratio, 
and P^[-] denotes a Fade approximant of order m in 
the numerator and 71 in the denominator. The 4 and 
5 PN coefficients, 05 and ag, are fit to numerical BBH 
waveforms. The values that give the optimal fit form a 
degenerate curve in the as-ag parameter space, and the 
specific values chosen here are (05, ag) = (0,-20). The 
second potential is rewritten as 

D{r) = B{r)A{r), (B3) 

and has been calculated to 2PN order 



D{u) = PI^[1 - Qvu^ + 2(3j/ - 2Q)vu% 



(B4) 



The motion of the EOB particle of mass ^ is deter- 
mined by the Hamiltonian 



V 



(B5) 



where 



off 
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Appendix B: Prescription for calculating EOB 
waveforms 



is the effective Hamiltonian. The equations of motion 
given this conservative Hamiltonian H and a dissipative 



In this appendix we compile the necessary ingredients 
needed to produce nonspinning BBH waveforms using 
the effective one body (EOB) formalism first introduced 
in Ref. The version used here is exactly that of 

Ref. [41] , and is described in more detail in a review [40] . 
The only ingredients not listed here are terms for the re- 
sumed waveform in Ref. 67 and coefficients to determine 
the ringdown waveform found in Ref. [68] . 



* The expressions below will be written exclusively in terms of 
rescaled dimensionless quantities. The coordinates (T, i?, 0) and 
conjugate momenta (Pi{, P^) have been rescaled to dimensionless 
coordinates {t,r,<j>) and momenta {pr,P4,) given by: t = T/M 
and r = R/M for the coordinates, and pr = Pn/ fi, = P^/fiM 
for the conjugate momenta. Other quantities are then rescaled in 
the following way: ui = MQ = Md<j)/dT is the angular velocity, 
D = D/M is the distance to the source, H = H/fi and H^ff = 
Hgf[/ii are the Hamiltonian and effective Ifamiltonian, and j^^ = 
T^/li is the radiation reaction force. 
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radiation-reaction force Ti are 



dr 
dt 


dH 

dpr 




(B7) 
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(B8) 


dpr 
lit 


dH 

dr 




(B9) 


dptj, 
~dt 


dH 

^~d6 




(BIO) 



Here, ^ = because the EOB Hamiltonian does not 
have an expUcit 4> dependence. In addition, for cir- 
cularized binary inspiral the radial component of the 
radiation-reaction force J> is of higher post-Newtonian 
order than the tangential component, so it is set to zero. 

To increase resolution near the black hole, the radial 
coordinate can be rewritten in terms of a tortoise coor- 
dinate [59] defined by 



dr 



1/2 



(Bll) 



The new radial momentum is then p^. = {A/ B) ' pr- 
Using this definition, the effective Hamiltonian becomes 



H 



eff 



Pi 



pI + A(l/r) 1 + -f + 2u{A - iv) 



(B12) 

where the parts that are 4PN and higher are neglected. 
(The 4 and 5 PN terms are however accounted for in the 
free parameters 05 and ag which were fit to numerical 
waveforms) . The equations of motion become 



dr 

dt 
d(j) 
lit 

dpr, 

dt 

dp^ 

~dt 



A dH 

VD dpr, 

dH 



A dH 



D dr 
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(B14) 

(B15) 
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Radiation reaction 



For the radiation reaction term J-^, which is written 
in terms of the PN parameter x, we will need a way to 
write X in terms of the dynamical variables. The usual 
method is to use the Newtonion potential 1/r and veloc- 
ity squared (wr)^ as PN counting parameters and then 
rewrite them in terms of the gauge invariant angular 
velocity lo using the Kepler law uj'^r^ = 1 which holds 
in the Newtonian limit, and for circular orbits, in the 



Schwarzschild {1/ 0) limit. The Kepler relation can be 
extended to circular orbits in the EOB metric by defining 
a new radial parameter, — r^p^^^, where 



i^{r,p^) 



dA 
dr 



l + 2v\ , A{r) 1 + 1^ I - 1 



(B17) 

for which oj^r^ = 1 holds for all circular orbits. In addi- 
tion, for noncircular orbits (in particular for the plunge) , 
this relation also relaxes the quasicircular condition by 
not requiring that the Kepler relation hold. The specific 
choice of PN parameter used here is 



X = [Lor^y. 



(B18) 



See Ref. [7D] for an extensive discussion. 

The radiation reaction term JF^ used in Ref. [JT] takes 
the form of a summation over all multipoles 



8 I 
1=1 m=l 



(B19) 



Instead of the standard Taylor expanded version of 
which can be found in Ref. [71] , Ref. [HZ] decomposed the 
waveform into a product of terms: 

h22 - h^r'S,sT22e''^- /22 {x)f^^'' (B20) 
for £ = TO = 2, and 

= /i?^-*4ffr,„e'*^"pL(a;) (B21) 

for the other values of i and to. The leading Newtonian 
part is given in the usual form as a function of x 



.Newt _ „ ^ (,,\^{i+<i)/2vt-e,-m 



IT 

2'' 



(B22) 



where the coefficients nim and ci+f_{v) are defined by 
Eqs. (5-7) of Ref. |S7|, and the parity e is for I + m 
even and 1 for i + m odd. 

The PN terms in the resummation which had been 
written as functions of x in Ref. |67| are now written in 
terms of the dynamical variables. The effective source 
term S^s becomes [72] 



<S'cff(r, Pr,,Pcl,) 

The tail term is 

Te,n{r,Pr,,P,p) -- 



Hcs{r,pr,,P4,) e = 
^ e = l 



r(^+l 2ik) j^ 2ik In 2kr 



(B23) 



T{£+1) 



(B24) 

where k = mnH{r,pr,,p^)uj{r,pr,,P(j,), k = 
muj{r,pr_^,p,j)), and tq — 2. The phase of this tail 
term is corrected with a term of the form e*''*™. The 
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first ten are given in Eqs. (20-29) of Ref. [57]. The 

first one is 

S22 = Iv"' + ^^V' - 2Auf'\ (B25) 

where y = {i'H{r,pr,TP^)u){r,pr,,p^)Y^^ and y, which 
lias several possible forms, is chosen to be y = w^/'^ [72] , 



Finally, the remainder term of the resummation fim is 
expanded in powers of x. For € = m = 2 this is then 
re-summed with a Fade approximant 

f22{x)^PM^'''°\x)l (B26) 

where 



,Tavior. x , 55i/ - 86 2M7 v"^ - Abv - 42%^ , 

fnn (y.x) = IH xA X 

J22 \ ) ^42 1512 



227875zy2 4l7r2jy 34625i/ 856 , , , , 21428357\ 

\ eulerhiolx) H x 

99792 33264 96 3696 105 ^ ' 727650 / 
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/ 36808 , , , , 5391582359 \ , / 458816 , , , , 93684531406 \ . , 
+ (, W^"^^'^^"^^") ~ 198648450 ) " + (,l984^^"^(") ~ 893918025 ) " ' ^^'^^ 

and the eulerln,„(x) = 7E + ln2 + ^ Inx + lnm terms are treated as coefhcients when calculating the Fade approximant. 
For the other values of i and m, /^m is re-summed in the form = p\ra- The quantity pi^n is given in Eqs. (C1-C35) 
of Ref. [67]. P21 is for example 

(IZv 59\ /617!/2 109931^ 47009\ , /7613184941 107 , , . x\ 1 

= l+(,^-56j"+(v^- Mm - 56448j" + (,2607897600 " 105 '^"^'^^^"^^"V " 

/6313 , , , , 1168617473883 \ , 

eulerlni x) x^ . B28) 

V5880 911303737344 / ^ ' 



The final product in the resummation of /122 is a next 
to quasicircular (NQC) correction term that is used to 
correct the dynamics and waveform amplitude during the 
plunge 



/2^2*^^ (01,02) = 1 + 



(ru;)2 



(B29) 



The free parameters a\ and a-i are determined by the 
following conditions: (i) the time when the orbital fre- 
quency w is a maximum (the FOB merger time Im^ co- 
incides with the time when the amplitude |/i22| is a max- 
imum, and (ii) the value of the maximum amplitude is 
equal to a fitting function that was fit to several BBH 
simulations and is given by 

|/l22|max(i^) = \.hlhv(\ - 0.131(1 - 4v)) . (B30) 



3. Integrating the equations of motion 

The equations of motion are solved by starting with 
initial conditions {rg, 0o,Pr,OiP0o} and numerically inte- 
grating the equations of motion. In this paper we are 
interested in long, zero-eccentricity orbits. This can be 
achieved in the FOB framework by starting the integra- 
tion with large r, where radiation reaction effects are 
small, and using the quasicircular condition p^, = 0. 



Fq. (B15l then becomes 
dH 



dr 



{r,Pr, = O,P0) = 



and results in the condition 



du 



{v?A{u)) 



(B31) 



(B32) 



for p0. If this quasicircular initial condition is used for 
smaller r, the radiation reaction term is no longer neg- 
ligable, and this initial condition will result in eccentric 
orbits. If desired, one can use an initial condition that 
more accurately approximates a zero eccentricity inspiral 
such as post-circular or post-post-circular initial condi- 
tions with nonzero Pr, |73j . 



To numerically solve Fqs. (B13 B16), they must be 
written as a svstem of fi rst order equations. However, 



the term J-,p in Fq. (B16) contains the square of f from 
the NQC term /^'^ in /i22- Since /^^^ gives a small 



correction of order 10% during the plunge, the easiest 
method, and that used in Ref. [H], is iteration [75]: (i) 
First solve the system of equations with set to one. 



(ii) Use the solution of Fqs. (B13 B16 1 to evaluate f and 
the other quantities in /^^"^^ . (iii) Re-solve the equations 
of motion with the NQC coefficients no longer set to one. 
(iv) Repeat steps (ii) and (iii) until the solution converges 
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to the desired accuracy. In practice this iteration only 
needs to be done roughly 2-5 times. 



A second method is to directly rewrite Eq. (BI6I as a 



first order equation. This can be done by replacing f in 
the NQC term on the right hand side with an expression 
containing and then solving for p^. The equations of 



motion (B13 BI6I and the chain rule give 




A dH 



L + M + Np^ 



where 



L 



M 




N ^ 



D dr dpl^ 
A d^H 
D dpr, dp^ 
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(B36) 
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Plugging this expression into Eq. (B16) yields an equa- 



tion quadratic in which can be solved exactly if 
desired. To first order in the NQC correction term, 



Eq. (BI6I now becomes the first order equation 

^QC 



.Higher 



22 



1 



+ 2f^{L + M) 



dt 
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P (B39) 



=2 m=l 
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includes just the higher order terms (£,m) ^ (2,2), and 
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The QC in means that the NQC term ^^2^^ has 
been factored out. 

The solution to the equations of motion 
{r{t), <j){t),pr^ {t ),P(i,{t)} are then plugged back into 



Eqs. (B20 -B21) to give the waveform h 



inspiral 
hn 



it). 



4. Ringdown 



In the EOB formalism the ringdown waveform of the 
final Kerr black hole is smoothly matched onto the in- 



spiral waveform at the EOB merger time Im- The mass 
of the black hole remnant is given by the energy of the 
EOB particle at the merger time tM 



A/bh = f^H{tM) = My 1 + MHesitAi) - 1), (B41) 
and the Kerr parameter is given by the final angular mo- 
mentum of the EOB particle [71] 



vp^{t 



M) 



M, 



BH 



1 + 2v{H,s{tM) ~ 1) 



(B42) 



The ringdown waveform is given by the first five posi- 
tive quasinormal modes (QNM) for a black hole of mass 
Mbh and spin aen: 



(B43) 



?i=0 



where cr — a22n + i'-^22n is the nth complex £ ~ 771 ~ 2 
QNM frequency for a Kerr BH with mass Mbh and spin 
obHj and are complex constants that determine the 
magnitude and phase of each QNM. The amplitude of 
the negative frequency modes is small |69j . The first 
three QNMs have been tabulated in Ref. [68] , and fitting 
formuli are also provided. The QNM frequency uj22n can 
be approximated by 



MBH'^22n = A + /2(1 - aBn)'' 



(B44) 



and the inverse damping time a22n is given in terms of 
the quality factor approximated by 



1 W22r; 



2 a22ri 



qi + 92(1 - flBH)'^^- 



(B45) 



The coefficients for n = 0-2 can be found in tabic VIII 
of Ref. [68] . For n = 3-4, Q!22n and a;22n can be linearly 
extrapolated from the values for n — 1 and 2 as was done 
in Ref. [73. 

The constants are determined by requiring that 

the inspiral and ringdown waveforms be continuous on a 
"matching comb" centered on the EOB merger time Im- 
Specifically, at the times {tM — 2S,tM — S,tM,tM + S,tM + 
26} we require hi^^'^''''\t) = /i^^down^^^^ ^ g 

was chosen to be equal to 2.3Mbh/M. This gives 5 com- 
plex equations for the 5 unknown complex coefhcients 

r+ 

^22n- 

The full inspiral plus ringdown waveform is then given 



by 



h22it) 



^ringdown f 
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\t) t>tM 
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